\[p(h, R) = \frac{2\pi R h}{2\pi R^2 + 2\pi Rh}\]
\[p(\eta) = \frac{\eta}{1 + \eta}\]
\[ p(h, R) = \frac{h}{2(2R) + h}\] \[ p(\eta) = \frac{\eta}{4 + \eta} \]
\[ p(\eta) = \frac{\arctan (\eta/2)}{\arctan(\eta/2) + \arctan(2/\eta)} \]
\[ p(\eta) = \,? \]
library(tidyverse)
n <- 60
m <- 38
mu <- 2
lambda <- 4
set.seed(47)
y_mu <- rpois(m, lambda = mu)
y_lambda <- rpois(n - m, lambda = lambda)
y <- c(y_mu, y_lambda)
df <- data.frame(t = 1:n,
y = y,
rate = rep(c("mu", "lambda"), c(m, n - m)))
ggplot(df, aes(x = t, y = y)) +
geom_point() + geom_line()
ggplot(df, aes(x = y)) +
geom_bar()
ggplot(df, aes(x = y, fill = rate)) +
geom_bar(position = "dodge2")
alpha <- 10
beta <- 4
nu <- 8
phi <- 2
it <- 5000
post_samples <- matrix(rep(NA, it * 3), ncol = 3)
colnames(post_samples) <- c("mu", "lambda", "m")
m_i <- n/2 # initialize m
for (i in 1:it) {
mu_i <- rgamma(1, alpha + sum(y[1:m_i]), beta + m_i)
lambda_i <- rgamma(1, nu + sum(y[(m_i+1):n]), phi + (n - m_i))
joint_vec <- rep(NA, n - 1)
for (j in 1:(n - 1)) {
joint_vec[j] <- mu_i^(alpha + sum(y[1:j]) - 1) *
exp(-(beta + j) * mu_i) *
lambda_i^(nu + sum(y[(j+1):n]) - 1) *
exp(-(phi + n - j) * lambda_i)
}
p <- joint_vec/sum(joint_vec)
m_i <- sample(1:(n - 1), size = 1, prob = p)
post_samples[i, "mu"] <- mu_i
post_samples[i, "lambda"] <- lambda_i
post_samples[i, "m"] <- m_i
}
joint_posterior <- data.frame(post_samples)
ggplot(joint_posterior, aes(x = mu)) +
geom_density()
ggplot(joint_posterior, aes(x = lambda)) +
geom_density()
ggplot(joint_posterior, aes(x = m)) +
geom_bar()
df <- data.frame(rates = c(post_samples[, "mu"], post_samples[, "lambda"]),
parameter = factor(rep(c("mu", "lambda"), c(it, it))))
ggplot(df, aes(x = rates, fill = parameter)) +
geom_density()